# Create a vector of package names
packages <- c("dplyr", "ggplot2", "tidyr", "tidyverse", "data.table", "readxl", "foreign", "devtools", "fixest", "sandwich", "modelsummary", "kableExtra")

installed_packages <- packages %in% rownames(installed.packages())
if (any(installed_packages == FALSE)) {
  install.packages(packages[!installed_packages])
}

# Packages loading
invisible(lapply(packages, library, character.only = TRUE))
rm(packages, installed_packages)

options(modelsummary_format_numeric_latex = "plain")
f <- function(x) format(x, digits = 3, nsmall = 2, scientific = FALSE)

getwd()
setwd(dir = "/Users/jmc4qg/The Lab Dropbox/Jonathan Colmer/ShotSpotter_env/Journal_submissions/ReStat/Replication Folder/")

mydata <- read.dta("Analysis Data/NIBRS_analysis.dta")
mydata_LP <- filter(mydata, MP == 0)
mydata_MP <- filter(mydata, MP == 1)

D <- feols(homicide_pc ~ 1, data=mydata)
D_LP <- feols(homicide_pc ~ 1, data=mydata_LP)
D_MP <- feols(homicide_pc ~ 1, data=mydata_MP)
D_diff <- feols(homicide_pc~MP, data=mydata)
models <- list(D, D_MP, D_LP, D_diff)
modelsummary(models, vcov = ~STATE, output = "Figures and Tables/Table_1/dstats_homicide_pc.tex", fmt = f, stars = c('*' = 0.1, '**' = 0.05, '***' = 0.01))

D <- feols(gun_homicide_pc~1, data=mydata)
D_LP <- feols(gun_homicide_pc~1, data=mydata_LP)
D_MP <- feols(gun_homicide_pc~1, data=mydata_MP)
D_diff <- feols(gun_homicide_pc~MP, data=mydata)
models <- list(D, D_MP, D_LP, D_diff)
modelsummary(models, vcov = ~STATE, output = "Figures and Tables/Table_1/dstats_gun_homicide_pc.tex", fmt = f, stars = c('*' = 0.1, '**' = 0.05, '***' = 0.01))

D <- feols(ngun_homicide_pc~1, data=mydata)
D_LP <- feols(ngun_homicide_pc~1, data=mydata_LP)
D_MP <- feols(ngun_homicide_pc~1, data=mydata_MP)
D_diff <- feols(ngun_homicide_pc~MP, data=mydata)
models <- list(D, D_MP, D_LP, D_diff)
modelsummary(models, vcov = ~STATE, output = "Figures and Tables/Table_1/dstats_ngun_homicide_pc.tex", fmt = f, stars = c('*' = 0.1, '**' = 0.05, '***' = 0.01))

D <- feols(agg_assault_pc~1, data=mydata)
D_LP <- feols(agg_assault_pc~1, data=mydata_LP)
D_MP <- feols(agg_assault_pc~1, data=mydata_MP)
D_diff <- feols(agg_assault_pc~MP, data=mydata)
models <- list(D, D_MP, D_LP, D_diff)
modelsummary(models, vcov = ~STATE, output = "Figures and Tables/Table_1/dstats_agg_assault_pc.tex", fmt = f, stars = c('*' = 0.1, '**' = 0.05, '***' = 0.01))

D <- feols(population~1, data=mydata)
D_LP <- feols(population~1, data=mydata_LP)
D_MP <- feols(population~1, data=mydata_MP)
D_diff <- feols(population~MP, data=mydata)
models <- list(D, D_MP, D_LP, D_diff)
modelsummary(models, vcov = ~STATE, output = "Figures and Tables/Table_1/dstats_population.tex", fmt = f, stars = c('*' = 0.1, '**' = 0.05, '***' = 0.01))

D <- feols(tMean~1, data=mydata)
D_LP <- feols(tMean~1, data=mydata_LP)
D_MP <- feols(tMean~1, data=mydata_MP)
D_diff <- feols(tMean~MP, data=mydata)
models <- list(D, D_MP, D_LP, D_diff)
modelsummary(models, vcov = ~STATE, output = "Figures and Tables/Table_1/dstats_tMean.tex", fmt = f, stars = c('*' = 0.1, '**' = 0.05, '***' = 0.01))

D <- feols(prec~1, data=mydata)
D_MP <- feols(prec~1, data=mydata_MP)
D_LP <- feols(prec~1, data=mydata_LP)
D_diff <- feols(prec~MP, data=mydata)
models <- list(D, D_MP, D_LP, D_diff)
modelsummary(models, vcov = ~STATE, output = "Figures and Tables/Table_1/dstats_precip.tex", fmt = f, stars = c('*' = 0.1, '**' = 0.05, '***' = 0.01))

